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TECHNICAL  REPORT  SUMMARY 

The  objective  of  this  project  is  to  combine  a  number  of  recent 
advances  in  finite  element  theory  and  computer  technology  for  analyzing 
cavities  and  structures  in  rock.  This  computer  program  applies  to  general 
three-dimensional  structures,  considers  nonlinear  material  properties, 
time  dependent  properties,  gravity  loading  and  sequence  of  excavation  or 
const  ruct i on . 

The  final  report  for  this  project  is  in  three  volumes  as  follows: 

Volume  1  -  "Analytical  Modeling  of  Rock  -  Structure 
Interaction."  Final  Report,  April  1973. 

Volume  2  -  "Users  Guide  for  a  Computer  Program  for  Analytic 
Modeling  of  Rock  -  Structure  Interaction."  Final 
Report ,  Apr i 1  1973 . 

Volume  3  ■  "Computer  Program  for  Analytic  Modeling  of  Rock  - 
Structure  Interaction."  Final  Report,  April  1973. 

All  three  volumes  were  prepared  by  Agbabian  Associates,  under  contract 
h0220035  with  U.  S.  Bureau  of  Mines.  The  project  was  sponsored  by  ARPA  under 
ARPA  Order  Number  1579,  Amend.  No.  3,  Program  Code  2F10. 

During  this  contract,  work  has  been  aimed  at  producing  a  user- 
oriented  computer  program.  The  work  of  writing  the  program  was  divided  into 
three  areas: 

a.  Input 

b.  Execution  and  output 

c.  Material  properties 


i  i  i 
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The  Input  Section  automatically  generates  the  continuum  part  of  the  finite 
element  mesh,  including  joint  elements,  allows  the  user  to  add  other  elements 
(beam,  shell,  truss)  to  the  mesh,  plots  the  result,  reduces  the  bandwidth  and 
reads  loads,  material  properties,  and  other  quantities  necessary  to  the  cal¬ 
culation.  The  Execution  Section  forms  the  global  stiffness  matrix  and  solves 
equations  of  equilibrium  for  displacements  by  an  implicit  method.  The  material 
properties  are  represented  by  subroutines  within  the  Execution  Section,  which 
are  written  in  a  modular  form  so  that  if  the  general  equations  of  nonlinear 
elasticity,  viscoelasticity,  viscoplasticity,  or  plasticity  do  not  suit  a 
particular  problem  they  may  be  easily  modified. 

One  of  the  guidelines  for  this  project  was  to  consolidate  existing 
finite  element  technology  into  a  single,  general  purpose  computer  program. 
Accordingly,  the  program  uses  existing  finite  elements,  a  proven  form  of 
the  equation  of  equilibrium,  existing  material  property  descr  i  p  t  i  or.s  ,  and 
an  existing  bandwidth  reducer.  However,  a  small  amount  of  new  work  was  done. 

A  new  joint  element  was  developed  and  an  existing  concept  for  automatic  mesh 
generation  was  greatly  extended.  Also,  a  form  of  Choleski  decomposition  was 
modified  for  efficient  use  of  mu  1 1 i buf fer i ng ,  resulting  in  substantial 

improvement  In  efficiency  of  solving  equations  of  equilibrium  using  peripheral 
storage . 

During  this  study,  some  new  work  was  reported  by  other  ARPA/Bureau 
of  Mines  contractors  which  has  been  Incorporated  in  the  program.  Among  these 
are  some  creep  data  obtained  by  W.  A.  Wawersik  of  the  University  of  Utah  and 
strength/deformabi lity  data  for  faults  by  R.  E.  Goodman,  F.  E.  Heuze,  anu 
Y.  Ohnishi  of  the  University  of  California,  Berkeley. 

The  technical  work  reported  in  Volume  1  is  divided  into  two  main 
parts.  The  first  part  describes  the  range  of  possible  application  of  the 
computer  program  to  mining  engineering  problems  and  reports  an  extensive 
analytic  study  of  the  Caladay  Project  hoist  room  with  which  experimental 
data  are  compared.  The  second  part  discusses  various  aspects  of  the  computer 
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program  and  its  theoretical  basis.  Attention  is  given  to  the  processing  of 
input  data  and  to  options  available  to  the  user  for  mesh  generation,  sequential 
excavation  or  construction,  automatic  bandwidth  reduction  and  plotting  of  the 
mesh.  Example  problems  which  have  been  solved  during  checkout  of  the  program 
are  described. 


Volume  2  explains  the  overall  operation  of  the  computer  program, 
defines  how  input  should  be  submitted  to  the  program,  defines  the  material 
property  models  in  detail,  explains  the  data  stored  on  each  tape  and  gives 
flow  charts  of  some  subroutines.  This  input  definition  is  for  the  version 
of  the  program  which  operates  on  the  UN1VAC  1108.  This  program  is  the 
primary  content  of  Volume  3,  and  is  available  on  seven  track  magnetic  tape 
from  DDC-TC,  U.  S.  Department  of  Commerce,  Springfield,  Virginia  22151, 
telephone  AC  (703)  321-8517.  The  tape  is  unlabeled,  even  parity,  external 
BCD  and  is  written  at  556  BP  I .  The  tape  has  constant  record  sizes,  each 
record  being  1920  characters  long  (24-80  column  card  Images  per  record). 

An  end  of  file  mark  follows  the  last  record  on  the  tape. 

Arrangements  to  obtain  copies  of  Volumes  2  and  3  may  also  be  made 
through  Agbabian  Associates,  250  North  Nash  Street,  El  Segundo,  California 
902A5,  telephone  AC  (213)  61*0-0576. 

This  contract  was  monitored  by  Dr.  William  J.  Karwoski ,  Spokane 
Mining  Research  Center,  U.  S.  Bureau  of  Mines. 
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FOREWORD 

The  principal  contributors  to  the  computer  program  described  in 
this  report  are  listed  below.  All  are  members  of  the  technical  staff  of 
Agbabian  Associates. 

A.  K.  Bhaumik 
K.  P.  Chuang 

J.  M.  Clark 

K.  T.  Dill 
J.  Ghaboussi 


J .  Isenberg 
E.  M.  Raney 
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INTRODUCTION 

The  main  purpose  of  this  document  is  to  show  how  the  input  is 
prepared  for  the  present  computer  program.  The  sequence  of  user  operations, 
logical  units  required,  format  for  input  data  and  flow  charts  for  some  of 
the  frequently-used  subroutines  are  given.  There  is  also  a  special  section 
on  the  input  for  material  properties.  The  prospective  user  should  refer  to 
the  final  report  entitled,  "Analytic  Modeling  of  Rock  -  Structure  Interaction," 
prepared  under  ARPA  Order  No.  1579,  Amend.  3  and  monitored  by  the  Spokane 
Mining  Research  Center,  US  Bureau  of  Mines  under  Contract  H0220035. 

Dr.  William  J.  Karwoski  was  the  project  monitor. 

This  computer  program  considers  static,  three-dimensional,  nonlinear 
analyses  by  the  finite  element  method.  Gravity  loading  and  sequence  of 
excavation  or  construction  are  also  considered.  The  following  types  of 
elements  are  available: 

a.  Beam 

b.  Rod 

c.  Ax i symmetric  quad 

d.  Plane  quad 

e.  Three-dimensional  hexahedron 

f .  Thick  she  1 1 

g.  Plane  slip  or  joint 

h.  Axisymmetric  slip  or  joint 

All  elements  except  the  beam  and  thick  shell  may  have  nonlinear  material 
properties.  The  types  of  nonlinear  properties  which  may  be  assigned  are: 

a .  Variable  modu \ i 


b. 


Elastic  or  variable  modulus  with  ideal  plasticity.  Yield 
function  may  be  isotropic  or  anisotropic 
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c.  Nonlinear  elastic  with  ideally  plastic  fracture  criterion 
and  strain  hardening  cap 

d.  Anisotropic  variable  moduli  with  anisotropic  fracture 
cr i terion 

e.  Viscoelasticity  (Kelvin,  Maxwell,  three-parameter  fluid). 

f.  Viscoplasticity 

The  material  properties  are  submitted  as  coefficients  in  an  array  whose  length 
is  100  decimal  locations.  The  coefficients  are  substituted  into  the  appropriate 
equations  according  to  flags.  The  user  may  easily  change  the  form  of  any  of 
the  equations  by  modifying  the  material  property  subroutines,  which  can  be 
done  independently  of  the  main  line  program.  Small  strains  and  small 
displacements  are  assumed.  However,  the  program  is  written  such  that  it 
can  be  modified  to  account  for  geometrical  nonlinearity.  The  main  provision 
in  this  regard  is  that  the  updated  coordinates  which  are  needed  to  form  strain/ 
displacement  relationships  and  the  original  relationships  are  available  in  the 
commutation  section  of  the  program. 

The  computer  program  was  originally  intended  to  solve  problems  in 
rock  mechanics  and  mining  engineering.  Thus,  the  excavation  and  construction 
options  are  intended  to  aid  analysis  of  shafts  or  chambers  in  rock  having 
nonlinear  properties  from  the  start  of  construction  through  excavation  and 
installation  of  temporary  supports  to  the  installation  of  permanent  supports. 
Besides  the  problems  shown  in  the  final  report  cited  above,  the  program  has 
been  applied  to  the  pullout  of  a  cable  anchor  embedded  in  soil.  This  problem 
uses  axisymmetric  elements  with  nonlinear  properties  and  axisymmetric  slip 
elements  for  slip  and  debonding. 

The  program  uses  the  direct  stiffness  method  of  structural  analysis. 
Degrees  of  freedom  are  defined  at  nodal  points  which  are  at  corners  of  elements. 
Element  stiffnesses  are  deposited  in  a  global  stiffness  matrix.  Nodal  point 
forces,  which  include  external  loads  and  internal  resisting  forces  are  expressed 
in  a  global  load  vector.  A  form  of  Choleski  decomposition,  involving 
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triangular  1 zat ion  of  the  global  stiffness  matrix,  reduction  of  the  load  vector 
and  back  substitution  of  the  reduced  lead  vector  into  the  triangularized 
stiffness  matrix,  is  used  to  obtain  incremental  displacements.  The  program 
has  several  options  for  performing  these  operations  from  which  it  chooses 
depending  on  the  available  core.  According  to  one  option,  the  stiffness 
matrix  is  formed  block  by  block  (one  bandwidth  wide)  and  stored  on  peripheral 
units;  it  is  then  retrieved  block  by  block  and  processed.  Under  a  second 
option,  the  stiffness  matrix  and  load  vector  are  formed  and  reduced  in  a 
block-by-block  manner  and  written  on  peripheral  storage  in  a  reduced  form. 

Part  of  the  reduction  is  performed  in  local  double  precision. 


These  operations  require  many  accesses  to  peripheral  storage.  The 
number  of  accesses  and  hence  the  computer  time  required  for  large  problems 
would  become  prohibitive  on  small  or  medium  sized  computers.  In  the  present 
program,  a  mul t i buffer ing  technique  is  used  which  allows  data  to  be  transferred 
to  and  from  peripheral  storage  while  computations  are  being  performed  in  the 
main  core  area.  This  procedure  virtually  eliminates  wait  time  for  Input/Output 
operations. 


The  computer  program  is  divided  into  three  sections  which  are  called 
Input,  Execution,  and  Output.  The  responsibility  of  the  user  is  greatest  in 
the  Input  section,  where  the  first  step  is  to  submit  by  means  of  coordinates 
on  cards  a  key  diagram  which  is  the  basis  of  the  finite  element  mesh.  The 
mesh,  tributary  areas  and  load  coefficients  are  generated  automatically,  and 
the  mesh  is  plotted.  At  this  point  the  program  will,  on  request,  generate  a 
restart  tape  and  terminate  execution  so  the  user  can  examine  the  mesh  he  has 
prepared.  Upon  restarting,  the  user  may  add  elements  or  nodal  points,  and 
he  may  request  the  bandwidth  reducer  to  attempt  reduction  of  the  bandwidth. 

At  this  point  a  new  plot  of  the  mesh,  which  may  have  been  expanded  or  renum¬ 
bered  or  both,  may  be  requested.  The  last  group  of  data  in  the  Input  section 
includes  element  activity  (an  element  may  initially  be  inactive,  become  active, 
and  then  become  inactive  again)  as  a  function  of  load  step,  load  data  (pres¬ 
sures  in  global  directions)  and  material  property  data.  At  this  point  two 
tapes,  a  linkage  file,  and  an  element  data  file,  linking  Input  to  the  Execution 
section  are  generated. 
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The  Execution  section  is  initialized  using  the  two  link  tapes 
created  b>  the  input  section.  The  user  supplies  the  number  of  steps  between 
reformulation  of  the  global  stiffness  and  the  load  step  after  which  a  restart 
file  is  to  be  generated. 

To  continue  execution  of  a  restarted  problem,  the  previous  restart 
tape  replaces  the  linkage  file,  and  the  run  is  continued.  Note  that 
reformulation  and  restart  intervals  may  be  changed  at  this  time. 

The  output  is  stored  on  a  magnetic  tape  and  is  printed,  on 
separate  request  by  the  user,  by  executing  the  Output  section.  Output  may 
be  requested  after  a  restart  has  been  generated,  and  after  the  run  is 
completed. 


SECTION  1 

DEFINITION  OF  INPUT  DATA 

The  input  data  cards  for  the  computer  program 
pages. 


are  described  on 


the 
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No.  of 

Cards 

Data  and  Definitions 

1.0 

OPTION  SELECTION 

The  following  cards  control  the  sequence  of 
operations.  Each  option  card  is  followed  by  the 
data  required  for  execution,  which  is  described 
in  subsequent  sections.  Options  may  be  stacked, 
but  only  in  the  sequence  given  below. 

1 .  *STAR 

2.  *ADD  data 

3.  *EXEC  data  5  3 

A.  *0UTP 

5.  *EXIT 


This  is  a  legal  sequence  allowing 

1 .  Mesh  generation 

2.  Modify  generated  mesh 

3.  Execute,  reforming  every  fifth  step  and 
generate  restart  at  Step  3 

**•  Print  stresses,  strains,  and 
d i splacements 

5.  Termination 

A  second  run  of  the  sequence: 

1.  •'■•REST  5  (N) 

2.  *0UTP 

3.  *  E  X I T 


AA 
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(cont) 


No.  of 
Cards 


Data  and  Definitions 


Fortran 
Forma  t 


Allows  restart  with  reformulation  every  fifth  step 
(optional  restart  at  Step  N)  followed  by  output 
and  termination. 


ISTR,  IWITIN,  NSTP,  NRSTRT ,  TIMER 

ISTR  =  control  character  identifies  this 
card  as  an  option  selection 

IWITIN  =  Four-character  abbreviation  for  the 
option  desired 

MODI  =  MODIFY — alter  contents  of  IWORDS 
array  (see  Card  Group  2.0) 

STAR  =  START--call  automesh  to  generate 
the  mesh  (see  Card  Group  3.0) 

ADD  =  Call  automesh  to  add  to  an  exist¬ 
ing  mesh  (see  Card  Group  A.0) 

EXEC  =  EXECUTE--initiate  calculation  (no 
additional  cards  required) 

REST  =  RESTART--restar t  calculation  (no 
additional  cards  required) 

OUTP  =  0UTPUT--begin  output  (no  addi¬ 
tional  cards  required) 

EXIT  =  End  of  run 


(A  1  ,A*t,5X, 
21 10.F10.0) 


NSTP,  NRSTRT,  TIMER  are  used  only  by  option  EXEC 
and  REST  and  should  be  blank  for  all  others. 

NSTP  =  step  count  for  reformulation  of  the 

global  stiffness,  [k]  will  be  reformed 
at  every  NSTP^  cycle,  must  be  >  0. 


AA  FCr<M  110.  107:> 
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Fortran 
Forma  t 

i.i 

(cont) 

NRSTRT  =  Restart  flag.  A  restart  tape  will  be 

generated  after  step  NRSTRT  and  calcula¬ 
tion  will  be  terminated. 

When  restart  is  requested,  the  code 

prints  a  message  indicating  which  file 

contains  restart  data.  A  run  may  be 

continued  from  this  point  by  mounting 
this  restart  tape  as  the  LINK  tape  (see 
Appendix  C)  and  requesting  option  *REST. 

TIMER  =  Restart  interval  in  seconds.  A  restart 

tape  will  be  generated  after  TIMER 

seconds  of  execution,  and  the  run  will 

continue.  Restart  tape  units  are  switchec 

among  available  units.  A  restart  unit  is 

not  released  until  a  new  restart  tape  has 

been  wr i tten. 

l 

aa  row;  no.  io/i> 
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Card 
G  roup 

2.0 


2.  I 


No.  of 
Ca  rds 


variable 


2.2 


2.3 


2.3.1 


Data  and  Definitions 


INPUT  FOR  ALTERING  IWORDS  VIA  OPTION  "MODI".  A 
DESCRIPTION  OF  IWORDS  IS  GIVEN  IN  APPENDIX  A. 

JTYPE  ,  N 

JTYPE  =  Flag  for  type  of  alteration 

-  I  ,  value  to  be  altered  is  integer 
(see  2.2) 

=  H  ,  value  to  be  altered  is  Holleritl 
(see  2.3) 

=  E  ,  end  of  data  for  MODI 
N  =  Number  of  cards  to  be  read 
Each  JTYPE  card  must  be  followed  by  N  cards  con¬ 
taining  data  used  to  reset  entries  in  the  IWORDS 
array.  These  cards  are  defined  in  Groups  2.2  and 
2.3. 

INPUT  CARDS  FOR  INTEGER  DATA  (JTYPE=I) 

JWD ,  IWD 

Set  IWORDS  (JWD)  =  IWD 

INPUT  CARDS  FOR  HOLLERITH  DATA  (JTYPE=H) 

JWD,  I HOL 

IF  (JWD  /  43),  set  I  WO  RDS  (JWD)  =  I  HOL 

IF  JWD  =  A3,  read  a  new  header  card  into 
IWORDS  (A3)- IWORDS  (62).  Card 
2.3.1  must  immediately  follow  the 
card  with  JWD  =  A3. 

(IWD(I),  I  =  A3 ,62) 


A  A 
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(215) 


(15,  IX, 

AA) 


(20AA) 


FC(<M  NO.  1075 


i  liait i  - 


. 


R— 72 1 5-2-270 1 


No.  of 
Cords 


DEFINITION  OF  INPUT  FOR 
BUREAU  OF  MINES  CODE 


Goto  and  Def i ni t i ons 


AUTOMESH  INPUT  DATA--0PT I  ON  "START" 

(HED(I),  I  =  1,  20) 

HED  =  A  80-column  label  of  alphameric  characters 
used  to  identify  the  problem.  This  label 
will  appear  at  the  top  of  each  printed 
output  page 

NSPANY ,  NSPANX ,  NSPANZ,  NMIDPT,  NCOMFS,  NCOMNP, 

IC00R,  NUMMAT,  IPLOT,  NSLIP,  NPELAD ,  ITSHEL, 

IRENIJM,  JPLOT ,  I  QUAD 

NSPANY,  NSPANX,  NSPANZ  =  Number  of  "spans"  or 
"blocks"  in  the  x,  y,  and  2  directions, 
respectively,  in  the  "key  diagram,"  NSPANZ 
is  not  meaningful  in  a  two-dimensional 
problem,  and  tne  corresponding  field  must  be 
left  blank. 

NCOMFS  =  Number  of  faces  (3-D  case)  or  sides 
(2-D  case)  to  be  formed  by  the  face-  or 
side-joining  process  described  in 
card  group  3-11  (page  16). 

NCOMNP  >  Number  of  pairs  of  nodal  points  which 
will  be  joined  b>  the  face-  or  side-joiring 
process'  in  Card  Group  3.11  (page  16). 

IC00R  =  A  numeric  flag  specifying  the  units  of 
length  in  which  the  input  data  is  given. 

=  0  if  the  uni ts  are  i n  feet 
“  '  the  uni  ts  are  in  inches 


Fort  ran 
Forma  t 


(ZOAb) 


(1515) 
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|  Data  and  Definitions 

If  units  selected  are  feet,  the  code  converts 
all  coordinates  and  areas  into  -inches. 
Material  property  data  is  not  converted  and 
is  expected  to  be  given  in  inches.  All  final 
mesh  data  is  expressed  in  inches. 

NUMMAT  =  Number  of  different  materials  in  the 
finite  element  mesh. 

'PLOT  =  A  numeric  flag  specifying  if  the  gen¬ 
erated  mesh  is  to  be  plotted. 

=  0  if  the  mesh  is  not  to  be  plotted 
=  1  if  the  mesh  is  to  be  plotted 

MSL'P  =  A  numeric  flag  signifying  the  presence 
of  slip  elements 

=0  if  no  slip  elements  are  in  the  mesh 
-  '  if  there  are  slip  elements  in  the'mesh 

When  slip  elements  are  present  in  a  mesh,  two 
runs  must  be  performed.  The  actual  number  of 
slip  elements  is  specified  only  in  the  second 

part  of  the  run  as  explained  in  Card  Group 

4.1. 


JNPELAD 


;  A  numeric  flag  specifying  if  there  are 
nodal  points  and/or  elements  to  be  man- 

Ua,,y  added  in  the  second  part  of  the  run. 
0  if  there  are  no  nodal  points  and/or 
elements  to  be  manually  added 


For t  ran 
Forma  t 
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=  1  if  there  are  nodal  points  and/or 

elements  to  be  added.  The- actual  number 
of  nodal  points  and/or  elements  to  be 
added  are  specified  In  Card  Group  4.1. 

ITSHEL  =  A  numeric  flag  signifying  if  there  are 
shell  elements  involved  in  the  mesh. 

=  0  if  there  are  no  shell  elements 

=1  if  there  are  shell  elements 

NOTE:  If  NSLIP,  NPELAD ,  or  ITSHEL  are  different 

from  zero,  a  restart  tape  will  be  generated 
and  option  ADD  must  be  run  to  complete  the 
mesh. 

IRENUM  =  A  numeric  flag  specifying  if  the  final 
mesh,  including  all  the  manually  added 
nodal  points,  elements,  and  slip  elements 
is  to  be  renumbered  for  possible  smaller 
bandwidth. 

=  0  if  the  mesh  is  not  to  be  renumbered. 

>  0  ,  IRENUM  minutes  will  be  allowed  for 

the  minimi zation 

JPLOT  =  A  numeric  flag  that  specifies  if  the 

renumbered  mesh  is  to  be  plotted.  This 
is  the  final  mesh  which  will  be  used  in 
the  numerical  calculation. 

=  0  if  the  final  mesh  is  not  to  be  plotted 
=  1  if  the  final  mesh  is  to  be  plotted 


Fortran 

Format 
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3.2  I  QUAD  =  This  variable  is  significant  only  in  a 

two-dimensional  run  in  which  the  mesh  is 
automatically  generated.  It  specifies 
if  the  quadrilaterals  generated  are 
plane  strain,  plane  stress,  or 
axi symmetric. 

=  -3  or  left  blank.  Quadrilaterals  are 
plane  strain. 

=  3  if  the  quadrilaterals  are  plane  stress. 

=  h  If  the  quadrilaterals  are  ax i symmetr ic. 

3.3  As  req'd  X  (l,  J)  (or  R ( I , J )  if  axisymmetric  mesh)  in  a  2-D  (8FI0.0) 

case  and 

X  ( I ,  J ,  K)  in  a  3*D  case 

X  =  X  coordinates  of  the  key  diagram.  Eight 

coordinates  are  given  on  a  card,  for  as  many 
cards  as  are  necessary  to  define  all  the 
nodal  points  in  the  key  diagram.  They  are 
given  in  the  order:  from  the  top  to  the 
bottom,  from  the  left  to  :he  right,  and  in 
the  3~D  case,  from  the  front  to  the  rear  of 
the  key  diagram. 

3.^  As  req'd  Y  (l,  J)  (or  Z ( I  ,  J)  if  axisymmetric  mesh)  in  a  (8FI0.0) 

2-D  case  and 

Y  ( I ,  J ,  K)  in  a  3-D  case 

Y  =  Y  coordinates  of  the  key  diagram.  Given  in 
the  order  as  stated  in  X. 
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As  req'd  Z  (I,  J,  K)  for  a  j-D  case  only 

Z  =  Z  coordinates  in  the  key  diagram.  Given  in 
the  order  as  stated  in  X. 

NMIDPT  I,  J,  X  (l,J),  Y(l  ,J)  in  a  2-D  case  and 

I.  J.  K,  X(  I,  J,  K),  Y  (I,  J,  K),  Z  (I,  J,  K) 
in  a  3_D  case 

I,  J,  K  =  positions,  i.e.,  row,  column  and  plane 
(or  layer)  numbers  of  a  midpoint  in  the  2-D 
or  3~D  array  of  nodal  points  and  midpoints 
in  the  key  diagram.  As  an  example,  in  the 
2-D  key  diagram  shown  below,  Midpoints  a 
and  b  occupy  the  positions  (3,6)  and  (*4,3) 
respectively. 


(8F10.0) 


ROW  1 - 

ROW  2 - 

ROW  3 - 

ROW  *4 - 


ROW  5 - 


COLUMN  1 


•MIDPOINT 


MIDPOINT  b 


3  **  5  6  7 


A  LINE  ON  WHICH  SOME  MIDPOINTS 
MAY  BE  SPECIFIED,  ALTHOUGH  NO 
MIDPOINTS  ARE  SPECIFIED  ON  THIS 
LINE  IN  THIS  KEY  DIAGRAM 


(21  5 » 2  F 1 0 . 0  ’ 
(3I5.3F10.0: 
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x,  Y,  z  =  CARTESIAN  COORDINATES  OF  THE  MIDPOINT 


3.7  As  req'd  ((MAT  (I,  J),  ,  ,  ,,  NSPANy)>  ,  .  NspsNX)  ()6|5) 

a  2-D  case  and 

((MATf  I  ,  J,  k),  1-1,  NSPANY) ,  J  =  1,  NSPANX)  ,  (1615) 

K  =  1 ,  NSPANZ)  in  a  3“D  case 

MAT  =  Material  number  of  a  block  in  the  key 
diagram,  16  numbers  per  card  for  as  many 
cards  as  are  necessary  for  all  blocks.  The 
material  numbers  of  the  blocks  are  given  by 
following  the  arrangement  of  the  blocks  from 
the  top  to  the  bottom,  from  the  left  to  the 


15 
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3.7 

(cont) 

right,  and  then  from  the  front  to  the  rear 
in  the  key  diagram.  The  material  numbers 
must  always  start  with  1,  and  a  zero  implies 
a  void  or  a  cutout. 

3.8 

As  req’d 

(NDY  (1) ,  |  =  1 ,  NSPANY) 

NDY  (1)  =  Number  of  divisions  the  blocks  located 

(1615)  ! 

in  the  1 th  row  of  the  key  diagram  are  to  be 
subdivided  into  in  the  y  direction,  16 
numbers  per  card. 

3.9 

As  req 1  d 

(NDX  (1 ) ,  I  =  1 ,  NSPANX) 

NDX  (1)  =  Number  of  divisions  the  blocks  located 

(1615) 

in  the  1  th  column  of  the  key  diagram  are  to 
be  subdivided  into  in  the  x  direction,  16 
numbers  per  card. 

3.10 

As  req’d 

(NDZ(  1),  I  =  1,  NSPANZ)  for  a  3-D  case  only 

NDZ  (1)  =  Number  of  divisions  the  blocks  in  the 

(1615)  ■ 

Ith  block  layer  are  to  be  subdivided  into 
in  the  z  direction,  16  numbers  per  card.  | 

3.11 

NCOMFS 

Nl, 

N2,  Ml,  M2  for  a  2-D  case  and 

(415) 

Nl, 

N2,  N3,  Ml,  M2,  M3  for  a  3-D  case 

In  a 

2-D  case : 

Nl, 

N2  =  The  end  points  of  a  mesh  line,  or  a  seg- 

ment  of  a  mesh  line,  which  is  to  be  fused 
with  another  mesh  line  with  end  points  Ml 
and  M2  as  defined  next.  Nl,  N2,  Ml,  and  M2 
are  temporary  nodal  point  numbers  that  are 
used  only  in  this  card  group  to  define  nodal 
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points  to  be  joined  in  the  key  diagram.  The 
nodal  points  are  numbered  from  the  top  to 
the  bottom  and  from  the  left  to  the  right  as 
illustrated  in  the  example  given  below. 


Ml,  M2  =  The  end  points  of  a  second  mesh  line, 
or  segment  of  a  mesh  line,  which  is  to  be 
fused  with  the  first  one.  Naturally,  there 
must  be  the  same  number  of  nodal  points  on 
this  line  as  on  the  first  line  defined  by 


N1  and  N2. 


The  two  mesh  lines  are  joined  or  fused  in  such  a 
manner  that  Ml  joins  N1  and  M2  joins  N2 . 


As  an  example,  a  rectangular  mesh  with  a  crack  or 
a  slip  surface  as  shown  in  the  figure  below  is 
the  result  of  using  the  accompanying  key  diagram 
and  specifying  that  mesh  line  2-6  is  to  be  fused 
with  mesh  line  3~7.  Naturally,  Nodal  Points  2, 

6,  and  10  in  the  key  diagram  have  to  be  given 
the  same  coordinates  as  Nodal  Points  3,  7,  and  11, 

respectively,  in  order  to  get  the  result  indicated 
in  the  final  mesh. 
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SCHEMATIC  KEY  DIAGRAM 

5 


1 

1  < 

1 

6 

0 

0  ^ 

,7  , 

1 

h- - i 

1  ^ 

i - < 

>TEMPORARY  NODAL 
10X  POINT  NUMBERS 

MATE  RIAL  NUMBERS 


ACTUAL  KEY  DIAGRAM  WITH  THE  CURRENT 
NODAL  POINT  COORDINATES 


10^. 

I  ^>SAME  AS 

II  NODAL  POINT 
COORDINATES 


If  it  is  specified  that 

1.  mesh  line  2-6  is  to  be  joined  to  mesh  line  3“ 7 
and  that 

2.  the  number  of  divisions  in  the  vertical  and 
the  horizontal  direction  are  to  be  2-1-1  and 
2-2,  respectively,  then  the  following  mesh  is 
obta i ned . 
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7  FINAL  NODAL  POINT  NUMBERS 
/  AUTOMATICALLY  ASSIGNED  BY 
19  THE  CODE 


|L5  V  120 
16  M2 1 


FINAL  ELEMENT  NUMBERS 


Note  that  variable  NCOMNP,  the  number  of  pairs  of 
nodes  to  be  joined,  must  be  specified  if  face 
joining  is  to  occur.  |n  the  above  example, 

Line  2-6  includes  three  node  points  (Nodes  3,  7, 
and  11  in  the  final  mesh)  so  that  NCOMNP  =  ? 
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3.11 

(cont) 


Data  and  Definitions 

In  the  3-D  example  below,  NCOMNP  would  equal  the 
total  number  of  nodes  generated  on  face  N1-N2-N3 
(which  must,  of  course,  equal  those  on  M1-M2-M3) 


In  a  3-D  case: 

N1 ,  N2 ,  N3  =  Three  corners  of  a  four-sided  face 
which  is  to  be  fused  wi th  another  face  as 
defined  by  the  corner  points  Ml,  M2,  and  M3, 
as  explained  below.  N 1 ,  N2  ...  M3  are 
temporary  nodal  point  numbers  that  are  used 
only  in  this  card  group  to  define  nodal  point: 
to  be  joined  in  the  key  diagram.  The  nodal 
points  are  numbered  from  the  top  to  the 
bottom,  from  the  left  to  the  right,  and  from 
the  front  to  the  rear  of  the  key  diagram. 

The  three  points,  N1 ,  N2 ,  and  N3  must  be  so 
specified  that  N1  represents  the  corner 
formed  by  the  two  sides  N1-N2  and  N 1  - N 3 »  as 
shown  in  the  following  illustration. 


Fortran 

Format 


Ml,  M2,  M3  =  Three  corners  of  a  second  four-sided 
face  which  is  to  be  fused  with  the  first  face 
as  described  above.  The  fusion  takes  place 
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(I  SIDE ( I )  ,  I  =  1  ,  A)  ,  NN  for  a  2-D  case  and 


( I  SIDE (I ) ,  I  =  1,  6),  NN  for  a  3-D  case 


In  a  2-D  case: 


Fortran 
Forma  t 


in  such  a  manner  that  Ml,  M2,  and  M3  coincide 
with  N 1 ,  N2,  and  N3,  respectively,  as  illus¬ 
trated  in  the  above  figure. 


(All,  6x , 
15) 

(611,  AX, 
15) 


ISIDE(l)  =  A  number,  1  through  A,  representing  any 
one  of  the  four  sides  of  a  key  diagram  for 
which  boundary  codes  are  to  be  automatically 
generated  by  the  program.  The  top  side  is 
designated  1,  the  bottom  2,  the  left  3,  and 
the  right  A,  and  these  numbers  may  be  given  in 
any  order.  Thus,  2A00  has  the  same  effect  as 
0A02;  both  of  them  say  that  boundary  codes 
are  to  be  automatically  generated  for  the 
bottom  and  the  right  side  of  the  key  diagram 
only.  It  must  be  emphasized  that  boundary 
codes  can  be  automatically  generated  only  for 
those  cases  in  which  the  sides  designated 
conform  to  the  following  requirements.  For  a 
side  of  the  mesh  which  corresponds  to  either 
the  top  or  the  bottom  side  of  the  key  diagram, 
the  side  must  be  parallel  to  the  x-axis,  and 
for  a  side  of  the  mesh  which  corresponds  to 
either  the  left  or  the  right  side  of  the  key 
diagram,  the  side  must  be  parallel  to  the 
y-axis.  A  blank  card,  however,  has  to  be 
supplied,  even  when  no  such  automatic  genera¬ 


tion  is  desi red. 
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In  a  3-D  case: 


Fort  ran 
Format 


ISIDE(l)  =  A  number,  1  through  6,  representing  any 
one  of  the  six  side  faces  of  a  key  diagram 
for  which  boundary  codes  are  to  be  automati¬ 
cally  generated  by  the  program.  The  front 
face  is  designated  1,  the  rear,  2,  the  top  3, 
the  bottom  A,  the  left  5,  and  the  right  6. 

All  the  rest  of  the  comments  given  above  for 
the  2-D  case  also  apply  to  this  case.  As 
before,  a  blank  card  must  be  supplied  even 
when  no  such  generation  is  desired. 

Both  in  a  2-D  and  a  3"D  case: 


NN 


NN  =  Total  number  of  nodal  points  for  which 

boundary  codes  are  to  be  read  in  manually. 

Refer  to  the  explanation  given  in  the  next 
card  group. 

J,  (KODE  (I,  J),  1=1,6)  (|5>  ix, 

611 ) 

J  =  Final  nodal  point  number  of  a  nodal  point  whose 
boundary  codes  are  to  be  manually  read  in. 

KODE  =  Boundary  condition  codes  for  the  six  degrees 
of  freedom  associated  with  the  nodal  point. 

=  0  when  the  degree  of  freedom  is  not 
constrai  ned. 


=  1  when  the  degree  of  freedom  is  constrained. 

Thus,  a  card  which  reads  25  011101  says  that  at 
Nodal  Point  25  only  the  x-b i splacement  and  she 
y-rotation  are  left  unconstrained. 
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3.13 

1 

(ISIDE 

(1),  1 

1  =  1 ,  A) ,  NN  for  a  2-D  case  and 

(ISIDE 

(1),  1 

1  =  1,  6),  NN  for  a  3-D  case 

Fortran 

Format 


(All,  6X , 
15) 

(611,  AX, 
15) 


Numbers  specified  in  this  card  group  parallel  those 


given  in  Card  Group  3.12.  The  only  difference  is 
that  the  numeric  flags  given  in  ISIDE  array  here 
designate  the  sides  or  faces  for  which  the  nodal 
tributary  areas  and  load  coefficients  are  to  be 
calculated  automatically,  and  NN  specifies  the 
total  number  of  nodal  points  for  which  tributary 
areas  and  load  coefficients  are  to  be  read  in 


manually  from  the  next  card  group.  All  other 
comments  given  in  Card  Group  3*12  apply.  As 
before,  a  blank  card  must  be  supplied  even  when 
no  such  automatic  generation  is  desired. 


(The  code  cannot  generate  tributary  areas  for 
axisymmetric  elements.) 


3. 1  A 


NN 


J,  (TA  (I,  J),  |  =  1,  3) 

J  =  Final  nodal  point  number. 


(I5.3F10.0) 


TA  =  Three  tributary  areas  associated  with  the 
three  components  of  the  nodal  force  acting 
at  Nodal  Point  J.  They  should  be  in  square 
ft,  if  the  coordinates  of  the  nodal  points 
in  the  key  diagram  have  been  specified  in 
feet.  In  the  code,  the  x-component  of  the 
nodal  force  acting  at  Nodal  Point  J,  for 
example,  is  defined  as 

Px  •  TA  (I,  J)  •  C (I ,  j) 
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UtMN  T  I  ON  OF  INPUT  FOR 
BUREAU  OF  HI  NET  CODE 


Data  and  Definitions 

As  an  example,  if  I S IDE  in  Card  Group  3.13 
read  3100,  or  0301,  or  0031  (as  long  as  3 
precedes  1,  that  is),  one  would  be  specifying 
the  following  variation  in  the  C  coefficient 
by  providing  the  following  two  cards  in  this 
card  group. 

First  card:  1.5  2.5 

Second  card:  -1.5 


Fortran 

Format 


© 

©  © 
© 


In  case  a  blank  card  is  supplied  in  Card 
Group  3.I3  (implying  0000),  no  card  in  this 
card  group  is  expected  by  the  program. 

In  a  3-D  case: 

CC  (1)  CC  {k)  =  The  values  of  the  load  coeffi¬ 
cient  at  the  four  corners  of  the  corresponding 
side  face  specified  in  Card  Group  3. 1 3 . 

All  the  comments  given  above  for  the  2-D  case 
apply.  The  following  fi}ures  define  the 
points  of  definition  of  the  four  coefficients 
for  each  side  face. 
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NN 


J, 


Data  and  Definitions 
(C(l,  J),  1=1,3) 


J  =  Final  nodal  point  number 

C  =  Three  load  coefficients  associated  with  the 
three  components  of  the  nodal  force  acting 
at  Nodal  Point  J. 


1  (VIEW  (I),  1=1,6)  for  a  3“D  case  only 

VIEW  (1  '  3)  =  Direction  cosines  of  the  horizon 

which  extends  from  the  left  toward  the  right. 

VIEW  (A  '  6)  =  Direction  cosines  of  the  view 
vector  (line  of  sight). 

The  plot  controls  specified  in  this  card  group 
are  for  the  plot  of  the  mesh  automatically  gen¬ 
erated  by  the  program.  No  cards  are  expected  in 
this  and  the  following  two  card  groups  in  case 
I  PLOT  =  0. 

1  QSIZE,  HTNODE ,  HTELE ,  NCUT 

QSIZE  =  Width  in  inches  of  the  paper  to  be  used 
if  different  from  the  default  value  of 
10.6  in. 

HTNODE  =  Height  in  inches  of  the  nodal  point 

numbers  if  different  from  the  default  value 
of  0.07  in. 

HTELE  =  Height  in  inches  of  the  element  numbers 
if  different  from  the  default  value  of 
0.08  in. 

If  IPL0T  =  0  this  card  is  not  supplied. 
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(6F10.0) 


(3F10.0, 

110) 
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Data  and  Definitions 


1  NCUT  =  Number  of  sections  to  be  defined  in  the 

next  card  group  which  will  define  the  portion 
of  thj  mesh  to  be  plotted. 

As  has  been  noted,  this  card  group  is  provided 
only  when  I  PLOT  =  1  . 

3.17.2  NCUT  ( (XCUT  (I,  J),  I  =  1,  6),  J  =  1,  NCUT) 

XCUT  (1  3,  J)  =  Cartesian  coordinates  of  a  point 

on  the  Jth  cut  plane. 

XCUT  (A  6,  J)  =  Direction  cosines  of  a  normal 
to  the  cut  plane.  The  normal  must  be  speci¬ 
fied  in  such  a  manner  that  it  points  toward 
the  side  of  the  cut  plane  which  is  to  be 
plotted. 

This  card  group  is  provided  only  when  IPLOT  =  1. 
For  a  run  which  does  not  contain  any  nodal  points, 
elements  or  slip  elements  to  be  manually  added, 
the  next  card  group  to  be  supplied  is  Card  Group 

4.5. 

AUTOMESH  RESTART  DATA  (OPTION  "ADD") 

NSLIP,  NNPACD ,  NELADD,  NOBB,  NUMMAT 

NSLIP  =  Number  of  generated  elements  to  be 

assigned  as  slip  elements  in  Card  Group  A. 2. 

NNPADD  =  Total  number  of  nodal  points  to  be  adder. 

NELADD  =  Total  number  of  elements  to  be  added. 

NOBB  =  Number  of  bar  and  beam  elements  to  be  added 
(also  included  in  NELADD). 
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Data  and  Definitions 


NUMMAT  =  Number  of  material  property  sets  to  be 
suppl i ed 

(NESLIP(N),  SLIPDA(l.H),  (IIJJ(J,N),  J  =  1,2), 
(SLIPDA(J,N),  J  =  2,3),  N  =  1,  NSLIP)  for  a  2-D 
case  and 

(NESLIP(N),  SLIPDA(l.N),  (IIJJ(J,N),  J  =  1,4), 
(SLIPDA(J.N),  J  =  2,3),  N  =  1 ,  NSLIP) 

NESLIP(N)  =  Element  number  of  an  existing  element 
(automatically  generated  in  the  first  run) 
which  is  to  be  associated  with  a  slip  elemen 
S  L I  PDA (1 , N )  =  A  numeric  flag  designating  the  side 
(2-D)  or  face  (3_D)  which  constitutes  part  o 
the  slip  surface.  The  following  figures  de¬ 
fine  the  numeric  designation.  Which  corner 
nodal  point  is  called  l,J,K,  etc.,  can  be 
determined  by  referring  to  the  printout  of  t 
first  run. 
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4.2 

1 1 J J  =  The  noda 1  po i n 

numbers  of  the  nodal  points 

(cont) 

which  are  located 

on  the  other  side  of  the 

slip  surface,  directly  opposite  to  the  nodal 

points  constituting  the  corners  of  the  side 

or  side  face  designated  by  ISLIP. 

In  the  case  of  a  2 

-D  problem,  there  are  two 

such  points  for  each  slip  element,  and  in  the 

case  of  a  3"D  problem,  there  are  four  such 

points  for  each  element.  These  additional 

nodal  points  must 

be  specified  in  the  follow- 

i ng  order  of  the  pairing  nodal  points: 

2-D  case: 

ISLIP  =  1 

1  ,  J 

ISLIP  =  2 

J,  K 

ISLIP  -  3 

K,  L 

ISLIP  =  4 

L ,  1 

3~D  case: 

ISLIP  =  1 

1,  J,  K,  L 

ISLIP  =  2 

J,  K,  KK,  JJ 

ISLIP  =  3 

K,  L,  LL,  KK 

ISLIP  =  4 

L,  1,  II,  LL 

ISLIP  =  5 

II,  JJ,  KK,  LL 

ISLIP  =  6 

1  ,  J,  JJ,  II 

SL 1  PDA  ( 2  ,  N )  and  SUPDA  (3,N) 

=  The  material 

number  and  the  thickness, 

respect i ve 1 y , 

of  the  associated  slin 

element  located  on  the  slip  surface. 
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I  TYPE ,  (IX  (I) 


Data  and  Definitions 
1=1,8),  MAT 


I  TYPE  =  Numeric  designation  of  the  element  type. 

-3  =  Quadrilateral  under  plane  strain 
condition 

-1  =  Bar  under  plane  strain  condition 

1  =  Bar  under  plane  stress  condition 

2  =  Beam  element 

3  =  Quadrilateral  under  plane  stress 

condi tion 

=  Axi symmetrical  element 

5  =  Brick  element 

6  =  Thick  shell  element 

7  =  Quadrilateral-slip  element 

8  =  Brick-slip  element 

9  =  Axi symmetric  slip  element 

IX  =  Element  connectivity  array.  The  nodal  points 
must  be  given  in  the  following  sequence: 

Bars:  Either  end  first 

Beams:  The  two  ends  and  then  the 
reference  point 
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Data  and  Definitions 


Quadrilaterals:  I,  J,  K,  L  in  the  counter¬ 

clockwise  sense  as  shown  in  the  figure 
on  Page  3&« 


Bricks:  I,  J,  K,  L,  II,  JJ,  KK,  LL,  as  shown 
in  the  second  figure  on  Page  36.  Notice 


that  when  the  brick  is  viewed  in  such  a 


manner  that  the  front  face,  l-J-K-L,  is 
closer  to  the  viewer,  l-J-K-L,  follow 
the  clockwise  sense,  and  that  II  on  the 
rear  face  is  located  on  the  same  edge 
as  I  ,  etc. 


Thick  shells:  The  16  nodal  points  involved 
are  specified  in  the  sequence  pictured 
in  the  following  figure: 


Quad-slip  and  axisym-slip  elements:  l,J, 
K,  L  for  the  quadrilateral  In  the 
sequence  already  explained,  and  then 
the  two  additional  nodal  points  as 
explained  under  Card  Group  4.2, 


I  tern  I  I JJ . 
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Data  and  Definitions 

Brick-slip  elements:  Eight  nodes  of  the 

brick,  I ,  J,  K,  L,  II,  J J ,  KK ,  LL 
followed  by  the  four  additional  nodes 
as  explained  under  Card  Group  A. 2, 

I  tern  I  I JJ  . 

Notice  that  when  the  number  of  nodal  points 
involved  is  less  than  eight,  some  fields  on  the 
card  are  skipped  before  MAT  is  specified.  However, 
when  the  number  of  nodal  points  exceeds  eight,  as 
in  a  thick-shell  element,  the  extra  nodal  points 
are  punched  on  an  additional  card  in  the  format 
(815).  The  read  i ng  of  this  extra  card  is  con¬ 
trolled  by  the  value  of  1TYPE. 

MAT  =  Material  number  of  the  element 

(KODE  (I),  I  -  1,  6),  X,  Y,  Z,  (TA  (l),  I  =  1,  3), 
(COEF  (I),  1=1,3) 

These  are  data  for  the  addition  of  nodal  points. 

KODE  =  Six  single-digit  numeric  flags  specifying 
the  constrained  condition  on  the  six  degrees 
of  freedom  at  the  nodal  point.  0  means  uncon¬ 
strained  and  1  means  constrained.  The  six 
degrees  of  freedom  are  arranged  in  the  follow¬ 
ing  order: 

1  x-d i splacement 

2  y-d  i  splacement 

3  z-d i spl acement 

A  x-rotation  (rotation  about  the  x-axis! 
•'•'Not  implemented  _  _ 


Fortran 
Forma  t 


(611 ,  2X, 
9F8.0) 
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5  y-rotation 

6  z-rotation 


x,  y,  z  =  Cartesian  coordinates  of  the  nodal  point, 


TA  =  Tributary  areas  associated  with  the  three 
components  of  the  external  pressure  acting 
at  the  point. 


C0Ec  =  Values  of  the  C  coefficient  associated 


with  the  three  external  pressure  components. 
The  x  component  of  the  nodal  force  acting  at 
a  point,  for  example,  is  eventually  computed 


where  P  is  the  x  component  of  the  externa 
X  I 


pressure. 


NACTIV 


NACTIV  =  Number  of  elements  with  activity  supplied. 


*4.5.1  NACTIV  NEL,  ION,  I0FF 


(315) 


NEL  =  Element  number 


Step  number  at  which  element  becomes  active 


I0FF  =  Step  number  at  which  element  becomes 


i nact i ve 


IDLFLG,  NDLINC,  (DCOSGR  (l),  1=1,3) 

IDLFLG  =  1,  output  is  relative  to  dead  loading 
=  2,  displacements  include  dead  loads 
NDLINC  =  Number  of  increments  in  which  to  apply 
dead  loads 

=  0  do  not  apply  dead  load 


(2  I  5  *  3F 1 0 . 


3*4 
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*4.6 

(cont) 

*4.7 

1 

*4.7.  1 

NLSTEP 

*4.8 

1 

*4.8. 1 

1 

*4.8.2 

NCUT 

*4.9 

*4.9.1 

1 

Data  and  Definitions 

DCOSGR  =  Coefficients  of  gravity  load  in  X,  Y, 
and  Z  d i rect i ons 

Dead-load  components  in  the  X,  Y,  and  Z  direction 
are  computed  via 


DC0SGRx  •' 

•  P  5 

vo  1  ume 

DC0SGRy  • 

volume 

DC0SGRz  - 

P  : 

vo  1  ume 

NLSTEP 

NLSTEP  =  Maximum  number  of  live-load  steps 
TIME,  PX,  PY,  PZ 

Time  =  Time  associated  with  load  step 

PX,  PY,  PZ  =  Pressure  loads  in  X,  Y,  and  Z 
d i rect i ons 

The  following  three  groups  of  cards  are  similar  to 
Card  Group  3-17-  They  provide  controls  for  the 
plotting  of  the  final  mesh,  after  any  nodal  point 
and/or  element  renumbering  has  been  performed. 

(VIEW  (l),  1=1,6)  for  a  3“D  case  only 
QSIZE,  HTNODE ,  HTELE ,  NCUT 

( (XCUT  ( I  ,  J) ,  I  =  1,6),  J  =  1,  NCUT) 

Material  Proper t i es--NUMMAT  sets  (see  Appendix  B) 
B,  G,  RHO 

B  =  Bulk  modulus 

G  =  Shear  modulus 

RHO  =  Dens i ty 


A  A 
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DEFINITION  OF  INPUT  FOR 
BUREAU  OF  HINES  CODE 


1 

m 

Data  and  Definitions 

Fortran 

Format 

4.9. 2 

i 

NISO,  NPL ,  NVISC,  ITER,  NSPLIT,  NITER 

(6110)  | 

NISO  =  Isotropic  flag 

NPL  =  Plasticity  flag 

NVISC  =  Viscoelastic  flag 

ITER  =  Iteration  flag 

NSPLIT  =  Number  of  splits 

NITER  =  Number  of  iterations 

4.9.3 

i 

DELTA,  CONV ,  NVM 

DELTA  =  Split  criterion 

CONV  =  Convergence  criterion 

NVM  =  Variable  modulus  flag 

(2F10.0.I 10 

4.9.4 

3 

(GCOEF  (1),  1  -  1,  20) 

(3 (8F10.0) ) 

Shear  modulus  coefficients 

4.9.5 

9 

(BCOEF  (1),  I  =  1,  40),  (YCOEF  (l),  1  =  1 

,  25) 

(7 (8F1 0. 0) , 
3F10.0.21 10, 
3F10.0/F10. 

BCOEF  =  Bulk  modulus  coefficients 

YCOEF  =  Inelasticity  coefficients 

4.10 

NOBB 

Additional  data  for  bars  and  beams 

( (ADBB  (1,  J),  1  =  1,  8),  J  =  1,  NOBB) 

(6F10.0, 

2110) 

ADBB  (1  8,  J)  =  The  additional  element 

data 

required  for  the  Jth  bar  or  beam  element. 

Cards  in  this  card  group  must  be  so 

arranged 

that  each  card  comes  up  in  the  same 

sequence 

as  the  corresponding  element  card  is 
tered  in  Card  Group  4.3. 

encoun- 
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Data  and  Definitions 

In  the  case  of  a  bar  element,  ADBB  (1,  J)  =  cross- 
sectional  area  of  the  bar. 

ADBB  (2,  J)  =  Initial  stress  in  the  bar. 

In  the  case  of  a  beam  element, 

ADBB  (1,  J)  =  A^ ,  the  cross-sectional  area  of  the 
beam. 

ADBB  (2,  J)  =  A^  the  shear  area  with  respect  to 
the  local  y-axis. 

ADBB  (3,  J)  =  A^,  the  shear  area  with  respect  to 
the  local  z-axis. 

ADBB  (A,  J)  =  |^,  the  second  moment  of  the  cross- 
sectional  area  with  respect  to  the  local 
y-axi s. 

ADBB  (5,  J)  =  1^,  the  moment  of  the  cross-sectional 
areas  with  respect  to  the  local  z-axis. 

ADBB  (6,  J)  =  I  or  J,  the  polar  moment  of  inertia 
xy 

of  the  cross-section.  If  this  field  is 
left  blank  J  will  be  computed  as 

J=l2+I3- 

ADBB  (7,  J)  and  ADBB  (8,  J)  =  The  two  end-release 
codes  for  the  beam. 


Fortran 

Format 
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APPENDIX  A 


CONTENTS  OF  IWD  COMMON  BLOCK 

Common  block  IWD  contains  constants  and  variables  used  throughout  the  pro¬ 
gram.  The  user  may  reset  the  default  values  listed  below  by  involving  option 
MODI"  (see  Card  Group  2.0  in  the  data  definition). 


Locat i on 

Name 

108 

Defaul t 

Value 

Meaning 

1 

IWCD 

5 

Card  reader  unit  number 

2 

IWPRT 

6 

Printer  unit  number 

3 

1 WPUNC 

-3 

Punch  uni t  number 

b 

IWC0MP 

1 108 

Computer  code 

5 

IWSRUN 

0 

Debug  flag 

6 

IWC0RE 

27000 

Length  of  BMCALC  common 

7 

IWDE0F 

4HDE0F 

End  of  file  word 

8 

IWBLNK 

Hollerith  code  for  blank 

9 

1  WIT  1  N 

(not  set) 

Current  itinerary  option 

10 

IWRUN(I) 

4HM0DI 

Allowable  itinerary  options 

11 

1 WRUN  (2) 

4HSTAR 

12 

1 WRUN  (3) 

3HADD 

13 

1  WRUN  C2*) 

4HEXEC 

lb 

1 WRUN (5) 

4HREST 

15 

1 WRUN (6) 

4H0UTP 

16 

1  WRUN  (7) 

4HEXIT 

17-19 

(Not  used) 

20 

IWNRUN 

10 

Length  of  IWRUN 

21 

IWTAPE 

8 

Number  of  tape  units 

22 

IWD RUM 

9 

Number  of  drum  units 

23 

IWLUN(I) 

9 

Tape  or  drum  logical  unit 
numbers 

2b 

IWLUN  (2) 

11 

25 

IWLUN  (3) 

13 

26 

IWLUN  (b) 

7 

27 

IWLUN  (5) 

12 

28 

IWLUN  (6) 

14 

29 

IWLUN (7) 

8 

30 

IWLUN  (8) 

10 

38 
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Location 

Name 

1108 

Defaul t 

V  a  1  ue 

Meaning 

31 

1 WLUN (9) 

21 

(Tape  or  drum  logical  unit 

32 

IWLUN(io) 

22 

number) 

33 

1 WLUN (11) 

38 

34 

1 WLUN (12) 

39 

35 

1 WLUN (13) 

4o 

36 

1 WLUN ( 14) 

41 

37 

1 WLUN ( 15) 

42 

38 

1 WLUN (16) 

hi 

39 

1 WLUN ( 1 7) 

44 

4o 

1 WLUN (18) 

0 

41 

1 WLUN (19) 

0 

42 

1 WLUN (20) 

0 

43-62 

IWHDR (l ) - 

Blanks 

Contain  run  title 

63 

IWHDR(20) 

IWNNP 

Not  set 

Number  of  node  points 

64 

IWNEL 

Not  set 

Number  of  elements 

65 

IWNMAT 

Not  set 

Number  of  materials 

66 

IWBAND 

Not  set 

Bandwi dth 

67 

IWNEQ 

Not  set 

Number  of  equations  (DOF's) 

68 

IWNLDB 

Not  set 

Number  of  entries/load 

69 

IWNREC 

Not  set 

coefficient  block 

Number  of  load  coefficient 

70 

IWNCF 

3 

blocks 

Number  of  columns  in  load 

71 

IWNEQB 

Not  set 

vector 

Number  of  equations/block 

72 

IWNEU 

Not  set 

Size  of  element  data  block 

73 

IWERR 

0 

Error  flag 

Ik 

IWLPPG 

50 

Maximum  lines/page 

75 

IWKEN 

27000 

Length  of  ATOMSH  common 

76-IOO 

IWDUM 

Not  used 

101 

LMSTRT 

21 

Pointer  to  start  of  LM  array 

102 

NRSTRT 

Not  set 

Restart  step  number 

39 
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Location 

Name 

1108 

Defaul t 
Value 

Mean i ng 

103 

TIMER 

Not  set 

Restart  time  interval 

104 

NSTP 

99999 

Reformulation  interval 

105 

MAX ST P 

Not  set 

Maximum  number  of  load  steps 

106 

NLDI 

Not  set 

Number  of  dead  load  increments 

107 

NCOL 

4 

Number  of  columns  in  load 
vector 

4o 
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A 

APPENDIX  B 

DEFINITION  OF  MmTFRIAL  PROPERTY  COEFFICIENTS 

The  coefficients  which  define  the  stress/strain  properties  for  each 
material  (N)  are  stored  in  an  array  COEFF(l.N),  I  =  1,100.  The  COEFF  array  is 
divided  into  three  main  sections  as  follows: 

1.  COEFF (I , N ) ,  I  =  1,15  Set  linear  elastic,  isotropic  moduli 

select  options  for  material  models, 
select  coefficients  controlling  inte¬ 
gration  of  stress/strain  relations 

2.  COEFF (I ,N) ,  I  =  16-35  Coefficients  defining  variable  shear 

modul i 

C0EFF(l,N),  I  -•  36-75  Coefficients  defining  variable  bulk 

modu I i 

3.  COEFF(I.N),  I  =  76-100  Coefficients  defining  inelastic 

models  (plasticity,  viscoplasticity, 
viscoelast ici ty) 

Although  the  properties  are  initially  read  into  the  COEFF  array,  they  are 
subsequently  transferred  into  dimensionless  variables  and  one  dimensional 
working  arrays  whose  names  suggest  their  meaning  as  shown  in  the  input 
definition.  Thus : 

Shear  Modulus(G)  COEFF(l.N),  I  =  16,35)  =  GCOF(l),  I  =  1,20 

Bulk  Modu 1  us (B)  COEFF(I.N),  I  =  36,75)  =  BCOF(l),  I  =  1,k0 

Inelasticity  Coef f i c i en t s  (Y)  COEFF(I.N),  I  =  76,100)  =  YC0F(l),  I  =  1,25 

To  make  writing  formulas  easier  in  the  following  discussion,  the  notation  G, 

B,  and  Y  is  used  to  indicate  the  various  coefficients.  Thus 

Gj  to  G£0  -  Shear  modulus  coefficients 

Bj  to  B^q  =  Bulk  modulus  coefficients 
Yj  to  Y^g  =  Inelasticity  coefficients 

M 


A 


R-72 1 5-2-2701 


The  dimensions  of  the  material  properties  are  expressed  in  the 
inch-pound-second  system. 

ELASTIC  MODULI  AND  CONTROL  COEFFICIENTS  (COEFF  ( I  ,N) ,  I  =  1,15) 

C0EFF(1,N)  =  Elastic  bulk  modulus;  constant  value  to  be  used  for 
isotropic  material.  May  be  used  with  plastic 
viscopiastic  or  viscoelastic  material. 

Elastic  shear  modulus;  constant  value  to  be  used 
for  isotropic  material.  May  be  used  with  plastic, 
viscopiastic  or  viscoelastic  material. 

C0EFF(3,N)  =  Mass  density 

(RHO) 

COEFF^.N)  =  Anisotropy  flag 

(NISO) 

=  0  Isotropic 

=  1  Elastic,  ideally  plastic  with  anisotropic  yield 

function  suggested  by  Hill 

=  2  Anisotropic  moduli  with  anisotropic  fracture 

cr i ter  ion 

Inelasticity  flag 
0  Elastic 

1  I dea i i y  plastic 

2  Strain  hardening  (cap  model) 

Viscous  property  flag 
0  Inviscid 
1  Viscopiastic 


COEFF (5 ,N) 
(NPL) 


C0EFF(6,N) 

(NVISC) 


COEFF (2 ,N) 
(G) 


2  Viscoelastic 
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C0EFF(7,N) 

(ITER) 


COEFF (8 ,N) 
(NSPLIT) 


COEFF (9 ,N) 
(NUMITR) 


COEFF ( 1 0 ,N) 
(SPLIT) 


=  Iteration  flag.  The  integration  of  nonlinear 
stress/strain  relations  may  require  iteration  if 
tangent  moduli  are  implicit  functions  of  stress 
and  strain  components.  This  is  especially  true 
of  the  cap  model  and  sometimes  true  of  variable 
modu  1  i  model s . 

=  0  For  iteration 

^  0  No  iteration 

=  Strain  subdivision  flag  for  viscoplastic  materials. 
To  compute  stress  relaxation  correctly,  the  strain 
increment  is  subdivided  into  NSPLIT  smaller 
increments.  This  number,  which  may  need  to  be  as 
high  as  50,  applies  to  viscoplastic  materials 
only. 

“  Number  of  iterations.  This  flag  is  used  with 
ITER  =  0.  Iteration  proceeds  until  convergence 
is  reached  within  a  prescribed  tolerance  (see 
CONV  below)  or  until  the  number  of  iterations 
equals  NUMITR.  May  be  about  8  to  10. 

=  Strain  subdivisior  flag  for  materials  other  than 
viscoplastic.  Subdivision  may  occur  when  a  trial 
stress  increment  exceeds  a  predetermined  value  or 
when  the  trial  stress  state  should  be  exactly  on 
a  yield  surface  but  if  off  by  greater  than  a  pre¬ 
determined  value. 

DELTA  0.  for  strain  subdivision  whenever 

do..'  .  ,  ■  DELTA.  In  this  case,  d  ..  are 

1  i  j  tr i al  i  j 

divided  into  DELTA  sub i nc remen ts . 


^3 


- - . . . 


R-72 15-2-2701 

DELTA  <  0.  for  strain  subdivision  whenever 

r  ^°i  j^t  rial 
h  =  ,  v  -  1 .  >  CONV 

(Atrial 

where 

f(a..)  .  ,  is  the  value  of  vJi  on  the 

i  j  t  r  i  a  1  *  2 

current  yield  surface,  and  &  i  s  the 

trial  value  of  the  second  invariant  of  the 

stress  deviator. 


In  this  case,  the  number  of  strain  sub i ncrements 
i  s 


COEFF (11 ,N) 
(CONV) 


(da .  . 
•J 

) 

*C0NV/F 

max 

trial 

where 

(da .  . 
•J 

^max 

is  the  largest  of  the 

trial 

trial 

stress  increments. 

To  avoid  subdividing  strain  increments,  set 
DELTA  » 

i  ii'  ma v i 

trial 


=  Criterion  for  convergence  of  iteration  (see  ITER). 
Also  used  to  control  strain  subdivisions  (see 
NSPLIT).  Convergence  is  defined  if 


a-b |  <  C0NV*a 


where 


a 

b 


idonl  +  I d022 1  +  lda33 

jdan  |  +  |da22|  +  |da33 


prev i ous 
step 


current 

step 
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COEFF ( 1 2 ,N) 

(nvm) 


COEFF ( 1 3  ,N) , 


=  Variable  modulus  flag 

=  0  For  constant  modul i  (uses  B  and  G  above) 

=  1  For  variable  moduli 

COEFF ( 1 ^ , N) ,  COEFF ( 1 5 , N )  are  presently  unused. 


SHEAR  MODULI  C0EFF(16,N)  through  COEFF(35,N)  (GC)F(1)  through  GC) F (20) ) 


Case  1 

(NISO  =  0,  NVM  =  0) 


(Isotropic,  constant  moduli) 
GC0F(1)  through  GCO F ( 20 )  =  0 


Case  ^  (Isotropic,  variable  moduli) 

(NISO  =  0,  NVM  =  1)  Shear  moduli  are  defined  to  be  a  function 

of  the  elastic  component  of  volumetric 
strain  u. 

di-n9  (vJ  <  umax  -  0)  where  ^max  is  the  smallest  previous 
va 1 ue  of  p . 


G  =  G2  "  (G2  '  Gl^  exP  (-M/G^) 
Unloading/Reloadinq  (p  <  p  <  0) 


Tension  (p  >  0) 


For  Case  2 

COEFF (N ,16)  =  G1 
COEFF (N » 1 7)  »  G2 
C0EFF(N,18)  =  G3 
C0EFF(N,25)  =  G1q 
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Case  3  (Constant,  anisotropic  moduli, 

(NISO  =  1,  NVM  =  0)  may  be  used  with  anisotropic  yield  criterion 

of  Hill). 

coeff(n , 1 6)  =  G12 
C0EFF(N,17)  -  G13 
coeff(n,i8)  -  g23 

(Variable  anisotropic  moduli  to  be  used  with 
anisotropic  fracture  criterion  of  Jaeger. 

May  be  used  for  plane  geometry  only.  Data 
on  fracture  criteria  must  be  specified 
through  JCOF.  Fracture  occurs  when 


Case  k 

(NISO  =  2,  NVM  =  1) 


elastic  (f  >  0) 

C0EFF(N,l6)  -  G12  (elastic) 

C0EFF(N,17)  -  G13  (elastic) 

COEFF (N, 18)  =  G23  (elastic) 

inelastic  (f  <  0) 

COEFF (N, 1 9)  =  G12  (inelastic) 

COEFF (N,20)  =  G13  (inelastic) 

C0EFF(N,21 )  =  G23  (inelastic) 


k  6 
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For  Case  2 

COEFF(N,36)  =  B 
COEFF (N , 37)  =  B2 

COEFF(N,38)  =  (presently  not  used) 

COEFF (N ,39)  = 

COEFF  (N,i*0)  -  B5 

C0EFF(H,1»5)  -  B10 


A 


Case  3 


R-72 1 5-2-2701 

(Constant  anisotropic  moduli.  May  be 


,  NVM  =  0) 

used 

wi  th  ani 

COEFF(N,36) 

=  B1 

= 

E1 

COEFF (N  ,37) 

=  B2 

= 

V12 

COEFF(N,38) 

=  B3 

= 

v13 

COEFF(N,39) 

=  Bl* 

= 

B2 

coeff(n.^o) 

=  B5 

= 

V21 

C0EFF(N,4l ) 

=  B6 

= 

v23 

COEFF (N , ^2) 

=  B7 

= 

E3 

C0EFF(N,i»3) 

=  B8 

= 

v31 

COEFF (N,M) 

=  B9 

= 

V32 

Case  4 

(NISO  =2,  NVM  =  1) 


(Variable  anisotropic  moduli  to  be  used  with 
anisotropic  fracture  criterion  of  Jaeger, 
may  be  used  for  plane  geometry  only.  Data 
on  fracture  criteria  must  be  specified 
through  YCOF.  Fracture  occurs  when 


f  =  (o1  -o^)  - 


2a0  ~  2  al°3 


-VT 


<  0 


+  1, 


elastic  (f  >  0) 


COEFF(N,36) 
COEFF (N , 37) 
COEFF(N,38) 
coeff(n,39) 
COEFF (N ,^0) 
C0EFF(N,41) 


-  B  =  E, 


=  B 


=  B 


=  B 


=  B 


=  v 


12 


=  v 


13 


=  B,  =  E, 


=  v 


21 


=  v 


23 


I 
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COEFF (M ,42) 
COEFF(N,43) 
COEFF (N ,44) 


E 

3 


inelastic  (f  <  0) 


COEFF (N,46) 

= 

B11 

= 

E1 

COEFF (N ,47) 

= 

B12 

= 

V12 

COEFF (N,48) 

= 

B13 

= 

v13 

C0EFF(M9) 

= 

B14 

= 

E2 

COEFF (N, 50) 

= 

B15 

= 

V21 

COEFF (N ,51 ) 

= 

Bl6 

= 

v23 

COEFF ( N , 52 ) 

= 

B17 

ss 

E3 

COEFF(N,53) 

= 

B18 

= 

V31 

coeff(n,54) 

= 

B19 

= 

v32 

COEFF(N,55)  through  C0EFF(N,65)  not  used. 


COEFF (N,66)  through  C0EFF(N,74)  contain 
direction  cosines  for  anisotropic  materials  as 
fol lows : 

GLOBAL  AXIS 


X 

Y 

z 

to 

—  >■ 

X  Q- 
<  O 

CtL 

X' 

e31 

% 

co 

CD 

<  O 

Q-  CO 

Y' 

B32 

B35 

B38 

O  Z 

Z  < 

cc  u. 

Q-  O 

V 

B33 

vD 

co 

CO 

cn 

CO 

CO 

^9 
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COEFF (N, 66) 
COEFF(N,67) 
COEFF(N,68) 
COEFF (N , 69) 
COEFF (N ,70) 
C0EFF(N,71) 
COEFF (N , 72) 
COEFF (N , 73) 
COEFF(N,74) ' 
COEFF(N,75) 


31 


32 


33 

B34 

B35 

B36 

J37 

J38 


39 

^40  =  (CCW  angle  from 

major  principal  stress  axis 
to  bedding  planes  at  which 
cohesion  of  anisotropic  rock 
is  a  minimum) 


COEFFICIENTS  DEFINING  INELASTICITY 

COEFF(N,76)  through  COEFF(N.IOO)  (YC0F(1)  through  YC0F(25)) 


Case  1 
(NPL  =  1) 


Elastic,  ideally  plastic  yield  criterion. 
Two  functions  are  available,  depending  on 
value  of  I  TYPE  =  C0EFF(N,95). 


For  I  TYPE  /  2,  material  may  be  isotropic  or  anisotropic  (Hill's 
criterion)  depending  on  values  of  coefficients. 
Principal  directions  of  anisotropy  are  submitted 
through  direction  cosines  in  BCOF.  For  ITYPE  =  2, 
criterion  is  that  used  with  cap  model,  but  cap  is 
not  available  unless  NPL  =  2. 


(NPL  1,  ITYPE  t  2)  Poly nominal  yield  surface,  with  von  Mises 

limit  at  high  pressure,  in  terms  of  J* 

Vi? 


50 


R-7215-2-2701 


J’l  ■  Yl°n  +  V2°22  +  v3°33 


\^2~ '■  “  Jr  Y4  (°lf  °22^  +  Y5  *°22  "  °ll)  +  Y6  ^°33 


2  2  2  I 

+  y  a  +  y  n  +  y  a  , 
7  13  8  23  9  12 


for  J,  >  Y10 


f,  =  Y„  ♦  Y,/  +  V13  (Jp2*Yl4(Jt)5-V^»0 


for  Ji  ;  Yio 


f2  ■  yi6  -  ^  >  0. 


(NPL  =  1 ,  I  TYPE  =  2) 


This  perfectly  plastic  yield  surface  is 
the  same  as  that  used  with  the  cap.  Cap 
is  not  used.  Anisotropy  is  not  recommended 

JL  |  I 

J"  **  and  J2*  ■=  J2;  hence,  Y1  through 
Yg  -  0. 

The  perfectly  plastic  yield  criterion  is 
prescribed  by 


for  J1  >  Y1q 


yi i  (i  -  h-^-l2)  +  y,,-  >py>o 


for  J,  <  Y,a 
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A 


Case  2 

(NPL  =  2,  I  TYPE  =  2) 


Cap  model  for  rock. 

Perfectly  plastic  yield  criterion  is  the 
same  as  that  given  above  for  NPL  =  1 ,, 

ITYPE  =  2.  In  addition  a  cap  is  prescribed 
by 


f2  =  (J1  -  V)2  +  P2  (J2  -  Q)  =  0 


where 


V  =  L  +  P  x  (L)  X1  (L) 

Q  =  I X ( L ) ]  2  {1  +  P2  [X  '  (L)  ]  2) 


for  L  <  Y, 


vn  -  n  -  y^i2  *  v 
10 


x 1  (l)  = 


- —  (l  -  L/Y  ) 

Y1q  U  L/T10; 


for  L  >  Y, 


X(L)  =  Y„  +  Y12 
X'  (L)  =  0. 


L  =  Y. 


f  (,  . 

J  o  \  .nr  .  r  / 


V3!- f, 


P  =  Y, 
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For  Cases  1  and  2. 


C0EFF(N,76)  =  Y 

COEFF(N,77)  =  Y2 
COEFF(N,78)  =  Y3 
C0EFF(N,79)  =  Y^ 

coeff(n,8o)  =  Y5 

COEFF (N , 8 1 )  =  Y, 

6 

C0EFF(N,82)  =  Y? 

COEFF(N,83)  =  Yg 

COEFF(N,84)  =  Yg 

C0EFF(N,85)  =  Y^ 

C0EFF(N,86)  =  y 

COEFF (N ,87 )  =  Y12 

COEFF (N,88)  *  Y 

C0EFF(N,89)  =  Yu 

C0EFF(N,90)  =  Y, _ 

I  P 

C0EFF(N,91)  =  y]6 

C0EFF(N, 92)  through  C0EFF(N,9i»)  not  used. 
COEFF(N,95)  =  ITYPE 


C0EFF(N,96)  through  COEFF(N.IOO)  not  used. 


Case  3 

(NPL  =  1,  NVISC  =  1) 


(Viscoplasticity.  May  be  used  only  with 
isotropic  moduli  and  an  isotropic  perfectly 
plastic  yield  criterion  of  the  type  described 
in  Case  1  (ITYPE  2)  above. 
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Y17  *(V 


3f 

3o .  . 

|J 


for  Case  3 

COEFF(N,76)  through  C0EFF(N,8^)  =  1.0 

COEFF(N,85)  =  Y]0 
COEFF(N,86)  = 

coeff(n,87)  =  y12 

C0EFF(N,88)  =  Y13 

C0EFF(N  ,89)  =  Y^ 

C0EFF(N,90)  =  Y15 

COEFF (N , 91 )  =  Yl6 

coeff(n , 92)  =  y1? 

COEFF(N,93)  through  C0EFF(N,100)  are  not  used. 
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Case 

(NISO 


For 


For 


(Anisotropic  fracture  criterion  of  Jaeger 
=  2,  NVM  =  1)  which  may  be  used  with  anisotropic  variabi.. 

moduli.  May  be  used  with  plane  geometry 
only.  The  fracture  criterion  is 


(oi 


2ao  ~  2YH.°3  , 


The  coefficients  a  and  a,  are  functions  of  fracture 
o  1 

parameters  which  are  submitted  through  YCOF.  The  appro¬ 
priate  parameters  are  selected  by  referring  to  two  angles 
as  foil ows : 

BETA  (computed  internally)  =  CCW  angle  from 

direction  of  major  principal  stress  to  the 
direction  of  the  bedding  planes. 

PSI  =  BCOF (hO)  =  Orientation  of  BETA  at  which 

a  is  a  minimum.  Usually  assumed  to  be 
o 

30°  (expressed  in  radians). 

Principal  directions  of  anisotropy  are  submitted  through 
direction  cosines  in  BCOF. 


0  <  PSI  <  BETA 


ao  =  Y11  "  Y12  lcos  2  (PSI  "  BETA)J 


BETA  <  PSI  <  90° 


a 

o 


17 


-  Y 


18 


[cos  2  (PSI  -  BETA)] 


19 
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1.0 


R-72 1 5-2-2701 


COEFF (N ,76)  through  COEFF(N,84)  = 
COEFF (N , 8 5 )  is  not  used . 

COEFF (N ,86)  =  Y 

COEFF (N , 87 )  =  Y12 

COEFF (N ,88)  =  Y13 

COEFF  (N  ,88)  = 

COEFF (N , 89) ,  C0EFF(N,30)  are  not  used 
COEFF (N , 9 1 )  =  yi? 

COEFF (N , 92 )  =  Ylg 

COEFF (N ,93)  =  Y 

C0EFF(N,94)  through  COEFF(N.IOO)  are  no 


used 


Case  5 
(NVISC  «  2) 


(Viscoelasticity.  May  be  used  only  with 
isotropic  moduli.  Kelvin  and  Maxwell 
elements  are  available.  Creep  may  be 
specified  either  in  volumetric  strain  or 
in  shear  strain  or  in  both  simultaneously 
depending  on  the  values  of  NVOL  and  NDEV. 


Kelvin  Element 


£t  exp  (-8,41)  +  ct  -1  (I  -  exp  (. 


At) 


creep  in  volumetric  strain  (NVOL  =  3) 
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creep  in  shf*or  strain 


a 


1 


2 


(NDEV  =  3) 

G 

nshear 

_ 1_ 

nshear 


Y 


Y 


7 


8 


Maxwell  Element 


at+4t  =  Et  +  (a24t  +  a1)cW  ■  Vt 

creep  volumetric  strain  (NVOL  =  2) 

a,  =  i  -  1/COEFF (1 ,N) 

a2  '  "bulk  '  V3 

creep  in  shear  strain  (NDEV  =  2) 

a}  =  1/G  =  1/COEFF (2 , N) 

a2  nshear  ^6 


3-Parameter  Fluid 

c 

Et+At 


't  +  Vt  +  a2° 


+  a^o 
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creep  in  bulk  (NVOL  =  1) 


2bul  k 


1  “  exp  (- 


K 


2bul  k 


At) 


1 


Ibul  k 


At  -  .  exp  (.  - At)\ 

'  n2bulk  ' 


n1bulk  +  n2bulk 


3  n1 bul k  K 


n 


At  - 


2bu  1  k 
K 


(l  -  exp  ( - - -  At)) 

'  n2bulk  / 


K 


2bu  1  k 


n1  bul  k  11 


1 


2bul  k 


=  Y. 


creep  in  shear  (NDEV  =  1) 


„  '2  shear 

1  G 


1  -  exp  (  - 


2shear 


At) 


1 


Ishear 


At  - 


2shear 


-  exp  (- 


2shear 


At) 


n1  shear  +  n2shear 


nl  shear  ^ 


n 


2shear 
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-  Y. 

n2shear 

- —  =  y 

n1shear  ^ 


For  Case  5 

COEFF(N,76),  COEFF(N,77)  not  used 

COEFF (N , 78)  =  Y3 

COEFF (N, 79)  = 

COEFF (N ,80)  =  Y$ 

COEFF (N ,8l )  =  Y, 

6 

COEFF ( N , 82 )  =  Y7 

COEFF (N ,83)  =  Yg 

COEFF (N ,8^4)  through  COEFF(NJOO)  not  used. 
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APPENDIX  C 
FILE  USAGE 


Logical  Unit 

Number  Name  Contents  and  Description 

9  LiNK  1.  Link  tape  contains  data  on  the  problem  size, 

loading  and  material  properties. 

Record  1  -  I WHDR (20)  -  run  title 

(38  words)  IWNNP  -  nunber  of  nodes 

iWNEL  -  number  of  elements 

IWNMAT  -  number  of  materials 

iWBAND  -  bandwidth 

IWNEQ  -  number  of  equations 

MAXSTP  -  number  of  live  loads 

NLDI  -  number  of  dead  load 

i ncrements 

1RELT  -  number  output/dead  load 
flag 

NTYPE(IO)  -  number  of  elements  of 
each  type 

Record  2  -  ID  array  (relates  node  points  to 

(6*1WNNP)  degrees  of  freedom) 

Record  3  “  Load  coefficients 

(3* iWNNP) 

Record  k  -  Live  loads-Time,  P  ,  P  ,  P 

x'  y  z 

(4-MAXSTP) 

Record  5  “  Material  properties 

(l 00* IWNMAT) 

2.  On  restarted  runs,  the  previously  generated 
restart  tape  is  mounted  as  LINK,  and  all  other 
files  remain  as  before. 
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A 

Logical  Unit 
Number 


Name 

LEU 


LOUT 


Contents  and  Description 

Element  data  tape  contains  one  record  for  each 
element.  All  records  are  the  same  size,  which  is 
the  maximum  length  required  for  the  elements  which 
are  present  (see  appendix  D  for  complete  descrip¬ 
tion  of  each  element). 


Output  data  tape. 


Record  1  -  IWHDR(20)  -  run  title 

(22  words) 

1WNNP  -  number  of  nodes 

iWNEL  -  number  of  elements 

Record  2  -  ID  array  relating  nodes  to  degrees  of 

(6*IWNNP)  freedom 


Record  3 
(*»  words) 


IWBAND  -  bandwidth 
iWNEQ,  -  total  DOFs 
I WNEQ.B  -  DOFs/block 
NBLOCK  -  number  of  blocks 


Records  k  -  N  contain  output  data  for  e£ch  load 
step,  written  as  required.  Each  record  is  IWNEQB+2 
or  23  words  long,  whichever  Is  greater.  The  first 
words  of  a  record  contains  the  time  corresponding 
to  this  load  step  and  word  2  is  a  signal  whether 
the  record  is  displacement  data  (if  «=  0)  or  element 
data  (*  number  of  elements  in  the  block).  If  the 
record  is  an  element  record,  each  element  consists  of: 

1.  element  number 

2.  element  type 

3.  yield  factor 

4.  stress  (6  words) 

5.  strain  (12  words) 

Element  data  is  packed  21  words  at  a  time  into  an 
output  block. 
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A 


Logical  Unit 
Number 


Name  Contents  and  Description 

LRST1  I  These  files  are  used  alternately  (as  required  by 
LFKI  /  the  whims  of  the  code)  to  contain:  l)  the  global 
LFKO  I  stiffness,  written  as  NBL0CK+1  records  of  IWBAND* 

NEQB  words  or  2)  Restart  data.  If  restart  is  used, 
the  fi le  must  hold: 

1.  65  words  of  common  data 

2.  70  words  of  control  data 

3.  20  words  of  array  addresses 

4.  Live  load  and  material  properties 
(4*MAXSTP  and  100*  IWMAT  words) 

5.  Global  stiffness  ( I WBAND*NEQB)  j  NBLOCK 

6.  Load/displacement  ( I WNCF*NEQB) J  sets 

7.  Stress  history  (< 248  words/e lement ,  IWNEL 
records) 

8.  Load  coefficients  (3*  number  of  loaded  Nodes  +  1 


20 

LSTRI  ) 

Stress  history  data,  written  1  record  per  element 

21 

LSTR0 j 

(at  most  248  words/element) 

t38 

LFIN  ) 

Contain  the  force  vector,  total  displacement,  and 

T39 

LF0UT j 

incremental  displacement.  NBLOCK  +  1  records 
1WNCF*IWNEQB  long 

40 

LOAD 

Load  coefficients,  three  words  for  each  loaded  node, 

written  in  blocks  of  NLDB  nodes 

+41 

LTMP1  ) 

Contain  element  stiffness  data  which  overflow  global 

142 

LTMP2 j 

stiffness  blocks.  At  mc-t  325  words  per  record, 

<  NUMEL  records 

+43- 

LTP1  ) 

Used  by  DEC0MP  to  store  terms  which  overflow  to  stiff. 

+44 

LTP2  j 

2 

ness  blocks.  Requires  at  most  IWBAND  +  NEQ  words 

These  files  are  not  needed  for  a  core  contained  (NBLOCK  ■  1)  probli 
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APPENDIX  D 
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CONTENTS  OF  ELEMENT  TAPE 


For  al 1  elements 


Length,  words 


IB  (1) 

IB  (2) 

IB  (3) 

IB  (b) 

IB  (5) 

IB  (6) 

IB  (7) 

IB  (8) -  I B  (11)  = 

IB  (12) 

IB  ( 1 3) “ I B  (21)  = 


Element  number 
Element  type 

Step  element  becomes  active 
Step  element  goes  inactive 
Minimum  DOF  for  this  element 
Maximum  DOF  for  this  element 
Material  number 
Volume  of  element* 

(Not  used) 

Di  recti  on  cosines 


b 

1 

9 


This  part  of  the  data  is  constant  for  all  elements.  Format  for 
remaining  data  is  similar  for  each  element. 

The  LM  array  always  starts  at  IB  (LMSTRT  +  l) 

Dead  loads  always  start  at  IB  (ND  +  LMSTRT  +  1) 

B-matrix  always  starts  at  IB  (2  *  ND+LMSTRT  +  1) 

LMSTRT  is  presently  set  to  21. 


Following  the  B-matrix  are  data  for  slip  elements. 


*IB  (8) - 1 B  (11)  contain  partial  volumes  for  axisymmetric  elements.  IB  (8) 
contains  total  volume  for  all  others. 
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Element  Type  +1 

1-D 

Bar 

ND 

IB  (22) 

DOF 

X1 

23 

DOF 

X2 

24 

DOF 

Y1 

25 

DOF 

Y2 

26 

DOF 

Z1 

27 

DOF 

Z2 

IB  (28) 

X1 

dead 

load 

29 

X2 

dead 

load 

30 

Y1 

dead 

load 

31 

Y2 

dead 

load 

32 

Z1 

dead 

load 

33 

Z2 

dead 

load 

IB  (34)  -  (39)  [B]  matrix  (6x1)  6 
IB  (40  -  45)  [B]  matrix  at  center  of  element  6 
IB  (46)  initial  stress  i 
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Element  Type  2  Beam  ND  =  12 


Length 


IB  (22-33) 


IB  (34-45) 


IB  (**6-189) 


DOF's 

xr  V  zr  °x.  ’  °y . *  ez1 

^9*  Y 9  »  ^0*  By  >  >  B-7 

*2  Y2  Z2 

Dead  load  components  corresponding 
to  DOF's 

Force-displacement  transform 


IB  (190-267)  -  Element  stiffness  [k] 
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A 

Element  Type  +3,  4  2-D  Quad  or  Axisym  ND  =  8 

IB  (22-29)  DOF's  Length 

*1'  *2’  X3’  V  V  V  Y3’  Yl|  8 

IB  (30-37)  Dead  load  components  8 

IB  (38-109)  [B]  matrix  as  3  x  24  72 

IB  (110-127)  [B]  at  center  of  element  as  3  x  6  18 

Element  Type  5  3-D  Brick  ND  =  24 

IB  (22-45)  DOF's  24 

Xj  1=1,8 

Y{  1-1,8 

z,  1-1.8 

IB  (46-69)  Dead  load  components  24 

IB  (70-333)  [B]  matrix,  or  3  x  11  x  8  288 

IB  (334-366)  [ B ]  at  element  center  33 
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ND  =  48 
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Element  Type  6  Thick  Shell 


IB  (22-69) 

IB  (70-117) 
IB  (118-405) 


NDOF '  s 

(X.,  Y.,  Z.)  1=1,16 

Dead  load  components 
e-u  transform,  or  6  x  48 


Length 

48 

48 

288 


Each  thick  shell  element  record  is  followed  by  five  records 
containing  [k]  for  the  element: 

IB  (1)  Highest  DOF  for  this  group,  ND1  1 
IB  (2- (ND1  +  D)  DOF's  ND1 
IB  ( (ND 1  +  2)-325)  Rows  of  [k] 


The  first  four  records  each  contain  six  rows  of  [k] ,  with  the 
last  24  rows  in  record  5. 

Values  of  ND1  for  the  five  records  are 


record  ND1 

1  48 

2  42 

3  36 

4  30 

5  24 
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IB  (22-29) 


IB  (30-33) 

IB  (3*»-*»1) 

IB  (*»2-*«5) 

IB  (146-117) 

IB  (118-135) 

IB  (136) 

IB  (137) 

IB  (138-1*41) 


IB  (1*46) 


DOF's  for  parent  element 


X,,  X2,  X3,  XV  V,.  V2,  Y  ^ 


DOF's  on  adjoining  slip 

Dead  loads  on  parent  element 
Dead  loads  on  adjoining  slip 

t B ]  for  parent  element 

t B ]  at  center  of  parent 

Face  number  of  slip  element 

Material  number  of  slip 

Coordinates  of  slip  surface 


V  X2’  V  Y2 


Slip  thickness 
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E'ement  Types  7  and  9  2-D  or  Axisym  w/slip  ND  =  12 


Length 

8 


72 

18 

1 

1 

*4 
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Element  Type  8 

JL 

3-D  with  SI  ip"  Nd  »  36 

Length 

IB  (12-45) 

DOF's  for  parent 

24 

X 

n 

CD 

CO 

II 

>- 

CO 

II 

IB  (46-57) 

DOF's  for  slip  surface 

12 

IB  (58-81) 

Dead  local  components  for  parent 

24 

IB  (83-93) 

Dead  load  components  for  slip 

12 

IB  (94-357) 

B  matrix  for  parent 

288 

IB  (358-390) 

B  at  center  of  parent 

33 

IB  (391) 

Face  number 

1 

IB  (392) 

Material  number 

1 

IB  (393-404) 

Coord inates  of  slip 

12 

IB  (405) 

Thickness 

1 

*Not  implemented 
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APPENDIX  E 


FLOW  CHART  OF  SUBROUTINE  TDRUM — 
PROCESS  ELEMENT  STIFFNESS  WHICH  OVERFLOWED 
PREVIOUS  GLOBAL  STIFFNESS  BLOCKS 
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c  TDRIM  ADDS  THE  OV 
ERF-LOtf  ELEMENT  K--BLOCKS  INT 
0  lit  GLOBAL  CP  30 


c 

i 

STIFFNESS.  FNE 

CFSSARY, 

OVERFLOW  IS  WRITE 

N  TO  THE 

CP  10 

- 1 - - - 

|c 

CURRENT  TEfPORARV  1 

FILE  FOR  FOLLOWING  BLOCKS. 

, - 

CP  50 

~  i  1 

COMMON  /IWD/IWCO,IWPRT,IW 
PLMC.  CWCOrP.IWSRLW.IWCORE 
.IWDC0C,:iM3L^,IWITIN/IWR 
IN  (10)  •IW'RN.MAPE.IWDR 
LMi  .EWLLM(20)  .IWHDR  (20)  ,IW 
M^P .  IWNEL .  .IWMMAT  >  IWBAMD.  I 
WNEQ.imJDB.IWhl^C.IWNCF. 
IWMLQB .  IWNEU « TVEER .  IWLPPG 

.:ewken.:ewdlh(2S)  ,lmstrt,n 

RSTRT .  FINER ,  NSTP .  MAXSTP .  N 
LDLNCOL 
— —  i 

COMMON  /BftlTSCAMiLEU.L 
OUT » LRSTl .  LFKE < LEKO .  LSTRO 
<  LSTRI ,  LF-1N .  LFOUT .  LOAD .  LT 
FPi  <  iTMP2  <  LTP1 .  L1P2 « NBLOC 
K.f^lJ.NTTPuNOlir.IVOUT.IO 
uT.NSTR.NLDB.NRFC.ftEGB.N 
EQBI.KFRRM.KCFNG.FIRST.NS 

teplftml.nel.nbltyp.nelma 

'TiND.LM.IB.MDFBLL.lsDPH^. 

DFIl*  NDF-0T2 » IEU»  IELK.  IKC  i 
NWF » NWF1N » ISTRDM » ISTFOT .  I 
STKW .  IWRT .  IST1 ,  ITMPl .  ITM3 
2,iyRK,IFI,ISCAtNWA.JRESB 
LIWRFCBMXM(7) 


DIMENSION  ELK(1),EU<2(1), 
AK(1) 


EQUIVALENCE  (D.NJ) 


R-721 5-2-2701 
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APPENDIX  F 

FLOW  CHART  OF  SUBROUTINE  BSTIF 
ASSEMBLES  GLOBAL  STIFFNESS  MATRIX 
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ca  m/Ai  cmtc  fKn  m 


R-721 5-2-2701 


c 

BEAM  ELEJ-ENTS 

AT  310 

© 

x^0T0\ 

JccrcmoNAJ — I?  ?> 


/  c 


■ifealiiB 


DWBM.H 


/  C 


■i.'IM'lli 


BLK(IELK)=D 


lLN=LEU 


WRITE  (IwPRT»120)IE0F 


CONTINLt 
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APPENDIX  G 


FLOW  CHART  OF  SUBROUTINE  ELPL-- 
COMPUTES  GENERALIZED  STRESS/STRAIN  RELATIONS 
AND  COMPUTES  STRESSES 


84 


START 


C  UM 


CE  10 


C  ELPL  CONTROLS  THE  GRO 

IP  OF  SUBROUTITJES  WICH  PRE 
PARE  THE  C  MATRIX  CE  50 

C  COMMUTE  TTE  STRESSES. 

.  BYPRODUCTS  OF  TEE  S£  COTP 
LfTATIONS  INCLLDE  E  CE  60 


C  AND  PLASTIC  COtfWT 
S  OF  STRAIN.  AW  MAXIMLM  EL 
ASHC  STRAINS. .  EL  CE  70 


C  ALSO  CONTROLS  TEERA!! 
ON  FOR  THE  APPROPRIATE  C  AN 
D  STRESSES,  AW  CO  CE  80 


C  TIE  CRITERION  FOR  SIR 
AIN  SPLITTING 

CE  90 


HM 


CE  100 


COMMON  /ELPLS/B.6.RH0.NIS 
O.NPL.NVISC,  NITER,  NSTEP.N 
ItttTR.  DELTA,  CONV.NVM.LHAX 
,VMAX,  MUTTER,  QQ.NSIZ 


COMMON  /BO/IWCO,IWPRT,IW 
PlNC,IWCOMP,IWSRlN,IWCORE 
,MEOF,IWBLT»<,MnN,IWR 
LM(iC)  ,IWNRLN,lWTAPE,IUDR 
IJM,IWLUN(20)  ,IWKDR(20)  ,IW 
WP,IWI€L,IW1AT,MAW,I 
^Q,mi3B,iyNREC,IWNCF, 
WWW.imimiyLPPG 
,IWKEN,IUD(11(25)  ,LMSTRT,N 
RSTRT,nre,NSP,MAXSTP,N 
LDI.NCOL 


© 
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DIMENSION  C(G.G)  ,Ci(G,G) , 
C2(6,6) , CEL  (6,6) , STRESS  (1 
2) ,STR0LD(6) ,STRSIN(6) ,EP 
SEL  (6)  ,STRINC(6)  .DEPSEL  (6 
)  ,DEPSPL(6)  ,E(12i«F(6),A( 
6)  ,H(6)  ,GC0F(20)  .BCOFHO) 
«YC0F(25) ,DIRC0S(3,3) ,DLM 
MY(i) , STRAIN  (12) 


c 

DEFINITION  OF  IWICF 

S' 

CE  220 

- 1 — - 

c 

NSIZ-  NUPER  OF  STRAll 

N  COfPOWNTS 

CE  230 

C 

NESO(=0  FOR  ISOTROP 

Y)  , 

(=i  FOR  ANISOTROPY) 

CE  210 

C  NVISC  (=0  F  INVISCID 
) ,  (=1  F  VISCOPUSTTC) ,  (= 
2  F  VISCOELASTIC)  CE  250 


C  NPL  (=0  FOR  EUSTIC) 
(=1  FOR  INELASTIC)  (=2  FOR 
CAPPED  MODEL)  CE  260 


C  NVM  (=0  FOR  CONSTANT 
MODULI!) ,  (-1  FOR  VARIABLE 
MOOULU)  CE  270 


C  NITER  (=0  FOR  ITERATE 
ON) ,  (=1  FOR  NO  ITERATION) 
CE  280 


C  NKETR  -  NllflER  OF  IT 
ERATIONS 

CE  290 


C  NELTYP  (=1  FOR  i-D  ST 
RESS) ,  (=-i  FOR  i-D  STRAIN) 

,  (=2F0R  BEAMS)  CE  300 


C  f€LTTP  (=3  FOR  2-0  ST 
RESS),  (=-3  FOR  2-0  STRAIN) 

,  (=8  FOR  SHELLS)  CE  310 


:  f&TTP  (=9 

FOR  AXISYfflETRIC  SUP) 
CE  310 


ITERION 


C  INITIALIZE 


| ALFA^ACOS (DIRCOS (1 .1) ) 


FACTORS 


CONTINUE 
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TEST  1 

F  (NVN.EQ.O.AfC.NISO.EQ. 

0) 

A\  ✓ 

1*)  T0\ 

<TEST  1  V — i  w 

nXlTTTfWAi  — /in 
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UKTIHi 


U  1*11.61. HHTTR) 


to  TO  X 
IfCITTCNAL 


IKUCmCNAL 


IC2(IJ)~*S+iCl(IJ)  H(IJ 


CONTTNUE 


